Nanoscale capillary wetting studied with dissipative particle dynamics 
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We demonstrate that Multi-Body Dissipative Particle Dynamics (MDPD) can be used as an 
efficient computational tool for the investigation of nanoscale capillary impregnation of confined 
geometries. As an essential prerequisite, a novel model for a solid-liquid interface in the framework 
of MDPD is introduced, with tunable wetting behaviour and thermal roughening to reduce artificial 
density- and temperature oscillations. Within this model, the impregnation dynamics of a water-like 
fluid into a nanoscale slit pore has been studied. Despite the coarse graining implied with the model 
fluid, a sufficient amount of non-equilibrium averaging can be achieved allowing for the extraction 
of useful information even from transient simulations, such as the dynamic apparent contact angle. 
Although it is found to determine the capillary driving completely, it cannot be intepreted as a 
simple function of the capillary number. 
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Introduction. Over the last decade, continuous, mcso- 
scale particle simulation methods such as Dissipative 
Particle Dynamics (DPD) 0,0] and many variants thereof 
El Q El have received considerable attention. Originally 
invented to include hydrodynamic effects in meso scale 
simulations of simple and complex fluids Q, 0, it also 
has successfully been_employed for studying polymeric 

lipid membranes [lOj, and 



systems or melts 0, 
for colloidal suspensions 
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As a model for solvents, 



one of the most important recent developments is the in- 
troduction of cohesive properties [lj, [lj, [l5| , extending 
the simple quadratic dependance on density in the equa- 
tion of state (EOS) in early formulations 0]. The ap- 
proach of Warren [TJj leads to particularly stable liquid- 
vapor interfaces and could be of great value in studying 
free-surface fluid dynamics problems where thermal cap- 
illary fluctuations are important, such as the intriguing 
phenomena related to the nanoscale Ray leig h instabil- 
ity (e.g. the break-up of liquid nano jets |lq|). In most 
numerical studies involving fluid particle (FP) ^3 meth- 
ods, the investigations have largely been kept generic. 
Specific interactions, e.g., between solid- liquid interfaces, 
have only been accounted for in crude ways, since in 
the majority of cases, the FP-intcractions do not arise 
from a systematic coarse graining procedure starting at 
the atomistic scale. We nevertheless suggest that a FP 
method can still be suitable for studying the aforemen- 
tioned phenomena, provided that adhesive and interfacial 
properties are introduced carfully and with respect to the 
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intrinsic molecular characteristics of a given FP model. 
It is the purpose of this paper to outline how this could 
be achieved, and we demonstrate for the capillary filling 
of a narrow slit pore that detailed information from short 
term transient simulations can be extracted. 

General Theory. Although the simulation scheme for 
FP-mcthods is MD-like in character, it employs two fea- 
tures notably different from traditional MD-simulations: 
soft interaction potentials of a finite range r c , and a mo- 
mentum conserving thermostat contributing a major part 
to the viscosity of the liquid 0,^3- ^ n DPD, particles in- 
teract via pairwise central forces Fy = F^ + F^ + F^- . If 
Ti denotes the particle position, the conservative force in 
standard DPD is F^ = Aw c (rij)e.ij, where = — Tj, 



and &ij — Tij/rij. The weight function 
(1 — r/r c ) vanishes for an inter-particle dis- 
tance r larger than a cutoff radius r c . The random and 



w c (r) 



dissipative forces are F£ 



3 and F ij 



— 7W D (r ) (vy • ey)e,j, respectively, and act as a ther- 
mostat if the amplitudes q of the random variable £y and 
the viscous dissipation 7 satisfy a fluctuation-dissipation 
theorem: q 2 = 2~fk B T and wfj(r) 2 = m°(r). The usual 
choice for the weight functions is w r 



iu G . Then, the 
EOS becomes at most quadratic in the FP-density p, ex- 
cluding the existence of capillary surfaces. However, F^j- 
may be augmented by density dependent contributions 
|l3lll4 Il5j in order to achieve arbitrary EOS. This class 
of schemes is termed multi-body DPD (MDPD). Here, 
the approach of Warren |l5j | is pursued who has the re- 
pulsive part of the force depend on a weighted average of 
the particle density, while the attractive part is density 
independent: 



P« = AijW^(rij)eij + + ^Mr^ey , (1) 



TABLE I: Parameters used in the simulations. The surface 
tension has been obtained in a planar geometry, by integrating 
the lateral part of the pressure tensor in normal direction. The 
characteristic fluid parameters are given in model units [m.u.] 
and in real units (MKS), where N m — 3. Note, that the time 
step At for the numerical integration of the MDPD equations 
is 2 orders of magnitude larger than At in a comparable MD 
simulation. 



Parameter 


Symbol 


m.u. 


MKS 


Fluid particle density 


P 


6.00 


(1000kg/m :i 


Interaction range (attr.) 


r c 


1.0 


(0.82nm) 


Interaction range (rep.) 




0.75 




Amplitude of F R 


q 


6.00 




Compress. 


dp/ dp 


45 ± 2 




Surface Tension 


a 


7.51 ±0.04 


(0.138N/m) 


Viscosity 


V 


7.41 ± 0.06 


(0.0003Pas) 


MDPD time step 


At 


0.02 


76 fs 



"The viscosity was determined by a fit to a Poiseuille velocity 
profile for a pressure driven flow in a slit of width b. 



with an additional weight function w d (r) = (1 — r/r^) 2 . 
For a single particle species, Ay = A is negative; the re- 
pulsive part with B > acts at a slightly smaller radius 
of interaction r<j. Pi &t the location of particle i is the in- 
stantaneously weighted average p { = J2i^j w d(rij)- War- 
ren showed that this approach produces stable capillary 
surfaces. In addition, we introduce a species dependent 
force constant A s ;, defining the interaction between liq- 
uid(l) particles and those of the solid(s) wall. The param- 
eters we use for this model (A = —40.0, B — 25.0) vary 
slightly from the ones in pjj. The length scale in real 
space is derived by matching the compressibility dp/ dp 
of the real fluid (water) to the one of the model fluid 
over a particle renormalization factor (the coarse grain- 
ing level) N m , as outlined in [icj. Formally, N m is the 
number of actual molecules a single fluid particle is to 
represent. The unit of time is set equal to r c /vth, where 
Vth means the (real) thermal velocity of one fluid particle 
with mass N m times the molar mass of a water molecule. 
This is independent of the transport properties of the 
system, as opposed to other prescriptions Q- All rele- 
vant quantities of our MDPD liquid arc summarised in 
Table [I] We stress that the case N m = 1 would not cor- 
respond to an atomistic limit; the MDPD liquid always 
remains a molecular substitute fluid, the molecular char- 
acteristics of which must be compared with the (contin- 
uum) phenomena studied. An adequate wall model. With 
FP approaches, there are roughly two classes of strate- 
gies to construct solid walls: (i) use of virtual boundary 
planes an d 1 ) walls made of frozen or crystalline 
arrangements of fluid particles [13, ED- In (i), some 
rule must be given to reintroduce a particle to the liq- 
uid when it crosses the boundary, wheras in (ii) the wall 
acts directly over interparticle forces. Models for sharp, 
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FIG. 1: (color online) (a) Density distribution within a 
half-slit (full width 20r c ) for a mutual attraction strength, 
A s i = —40.0 for a rigid wall (red line) and thermally rough- 
ened wall (green line). For a rigid wall unphysical density os- 
cillations occur near the wall, (b) Variation of contact angle 
with varying A s i , as determined from a simulated optical mea- 
surement (•), or using Young's law: a sv — o s i + <?i v cos(#o), 
and obtaining <r s ; from integrating the lateral pressure tensor 
in planar geometry (■). 



rigid interfaces will usually produce order effects and thus 
density oscillations that on the targeted length scale are 
not expected. Also, inhomogenous temperature profiles 
may arise |22|, or other s pur ious effects (detachment of 
droplets under shear flow [21|). In this work, we shall use 
a thermally roughened wall and allow fluid particles to 
penetrate the solid-liquid interface slightly. To achieve 
thermodynamic compatibility, the solid phase is made of 
fluid particles at the same density as the fluid, result- 
ing in a strictly homogeneous temperature profile across 
the interface. The wall particles are pinned by harmonic 
forces to their initial position, and mutually interact via 
ordinary MDPD forces (eq. (1) with A ss = -40.00). An 
additional, soft harmonic force acting normal to the in- 
terface prevents single fluid particles from diffusing into 
the wall indcfinctcly. Fig. ^ illustrates density profiles 
and the variation of the static contact angle of this solid- 
liquid interface. In contrast to a conventional rigid wall 
density oscillations can be much reduced (compare red 
and green lines in Fig. la, corresponding to two differ- 
ent strength of the pinning force) . Note that the wetting 
behaviour of the wall can be tuned arbitrarily between 
hydrophilic and hydrophobic by diminishing the mutual 
interaction A s i between solid- and liquid particles. The 
interfacial energy can be measured with respect to a pla- 
nar interface, by integrating over the transverse part of 
the pressure tensor; the static contact angle 9o resulting 
from the Young-Laplace equation is in excellent agree- 
ment with an optical measurment (cmp. Fig. lb), where 
9q is determined from the meniscus shape of a liquid slab 
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trapped between parallel walls. Details about this model 
of a solid-liquid interface will be published elsewhere • 
Dynamic behaviour of the impregnation process. In 
order to demonstrate the applicability of our interface 
model to nanoscale impregnation, we study the dynam- 
ics of a water-like liquid in a slit pore of width b = 20 
nm ~ 24r c attached to a finite reservoir (see Fig. 2a for 
technical details). We consider the case where 9o = 0, 
(i.e. A s i = —40.0). Fig. 2b and c depict snapshots of 
the meniscus at two different times during impregnation. 
Partly due to pronounced capillary oscillations that are 
excited along with the release in free energy during im- 
pregnation, the meniscus deviates from a spherical shape. 
By ignoring contributions of the wedges, an apparent dy- 
namic contact angle 9&{t) can be extracted by a fit to a 
circle segment (red dashed lines in Fig. 2b, c). The evo- 
lution of cos(#d(t)) and the penetrated height z{t) over 
time (red solid lines in Fig. 2d, e) can be divided into 
two regimes. While for t < 0.35 ns the liquid motion is 
dominated by inertia and high velocities, resulting in a 
rapid increase of 8d, a slow decrease towards the static 
contact angle can be observed for later times. It is impor- 
tant to note that the spatial dimensions and time scales 
probed here (the hydrodynamic time scale Th and the 
relevant capillary oscillation period t osc ) are still well 
above the intrinsic molecular scales of the model fluid 
|26j . This suggests that phenomenologically, the impreg- 
nation dynamics should follow some continuum law (as 
it should for a fully atomistic description) balancing vis- 
cous, capillary and inertial forces with the momentum 
change d/dt(m(t)z(t) + am res (t)z res (t)) . Here the time 
dependent velocities (i(i), z res (t) and accelerated masses 
(m(t),m res (t)) of the liquid in the slit and a fraction a 
of the liquid in the reservoir have been introduced. For 
a slit geometry as given in Fig. the evolution of the 
height of the miniscus is governed by 

pz(c.\z + al ) + pc 1 z 2 = cos(#d) pr^i, (2) 

o tr 

with ci = (1 — ab/b res ), where b/b res denotes the ratio 
of the slit- to the reservoir width in a;-direction. The 
Laplace pressure on the right hand of Eq. 2 side is the 
driving force of the impregnation and consequently one 
expects the dynamics solely to be determined by the time 
dependent behaviour of the (apparent) dynamic contact 
angle 6>d. This is indeed the case: if the measured course 
of cos(#d(i)) (i'cd solid line in Fig. 2d) is inserted into 
Eq. [3 a numerical integration results in the dashed- 
dotted line in Fig. 2e, and is in excellent agreement with 
the numerically measured z(t), with fluctuations sup- 
pressed by viscosity and inertia. As we may speak of an 
average 6^, it would now be appropriate to make contact 
to standard continuum theories using constitutive laws of 
the form #d=#d(Ca), where here Ca = rjz/a denotes the 
capillary number. The straightforward insertion of em- 
pirical laws that have been obtained experimentally for a 
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FIG. 2: (color online) Impregnation of a slit pore. While 
the simulation cell is finite in z— direction, periodic boundary 
conditions are applied in x and y. The period in y-direction 
(slit depth) is set to 12r c . The simulation involved ~ 170, 000 
FP, a complete run requiring 36 hrs. on a single proces- 
sor workstation, (a) Initially, a spherical meniscus with the 
static contact angle 9o = 0° has been prepared. Panels b and 
c show snapshots of the dynamic contact angle at later times. 
Panel d compares the MDPD impregnation dynamics (solid 
line) with corresppnding predictions from various continuum 
models. The dash-dotted line results from an analogous inte- 
gration, with the empirical law of Ref. |2^l . Triangles (green) 
represent 008(0(1) (inset) and z(t), with 6d determined from 
a stationary MDPD plug flow experiment. Dashed (green): 
integration of Eq. [5] with driving Laplace-pressure derived 
from the measured dynamic contact angle 8d, see inset for 
cos(#d(i))- As a general reference, the Washburn solution |25ll . 
valid for asymptotically large times (Q& ~ do), is displayed. 



similar range of Ca as found here (Ca ~ — 0.13) fails in 
general (see the law of Joos and Bracke [24| as an exam- 
ple, Fig. 2c, dashed line). Analysing #d(Ca) in terms of 
theoretical models such as a chemical or a hydrodynamic 
model [2?! or combinations thereof [2^, however, is pos- 
sible but appears not to be very conclusive. Rather, it 
is instructive to extract a constitutive law #d(Ca) from a 
steady state MDPD simulation with stationary flow pro- 
file. For the relevant range of Ca, a sufficiently long fluid 
plug is driven at constant speed v through an infinite 
slit pore with the same width as in the impregnation 
study. The capillary oscillations are much reduced in 
this setup, and the meniscus profile and 9^ can extracted 
in a co-moving frame. If we then insert the resulting 
6>d(Ca(v = z(t))) into Eq. the deviation from the mea- 
sured impregnation is again significant (green triangles 
in Fig. 2d,e). 
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Discussion and conclusions. The apparent "failure" of 
the consistency check above is rooted in the implicit as- 
sumption that in general, the (apparent) dynamic con- 
tact angle is a function of Ca and possibly a small set 
of parameters: local dissipation mechanisms in a liquid 
wedge close to or at the triple phase contact line are as- 
sumed to govern its mobility [27j • However, the realistic 
numerical setup (afforded by the numerical benefits of 
the FP method) and the molecular nature of our solid- 
liquid model system include aspects that go beyond a 
simple picture: (i) as the simulation is full 3D, the con- 
tact line is not only allowed to undulate in x, but also in 
y— direction. This should affect its average mobility [jfij . 
(ii) More important, our results suggests that #d depends 
on the flow field as a whole, as the Poiscuillc profile of 
the stationary experiment and the approximate plug flow 
in the inertial phase of the impregnation process must be 
opposed |33| . A functional dependance of 9d on the flow 
field is discussed in |2flj| : a key feature is the dynamics of 
interface formation, allowing for nonequilibrium values 
and gradients in surface and interfacial tension. 

It must be emphasized that since the model fluid does 
not employ atomistic interactions (e.g., of the LJ-typc), 
the continuum aspects of the "true" impregnation dynam- 
ics are possibly amplified. Nevertheless, we would like 
to stress the usefulness of the FP approach as an effi- 
cient tool for the exploration of relevant processes close to 
and at the contact line. In combination with subsequent, 
careful large scale atomistic verfication, FP methods can 
provide important contributions to a unifying picture of 
the rich physics occuring in dynamic wetting phenomena. 
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